The GVT space-time method for the self-energy of large systems 



Martin M. Riegeri, L. Steinbeckt *, I. D. White§, H. N. Rojas§0, and R. W. GodbytS 

^ Cavendish Laboratory, University of Cambridge, Madingley Road, Cambridge, CBS OHE, UK 
''Department of Physics, University of York, Heslington, York YOl 5DD, UK 

(February 1, 2008) 



OO 

0\ 



o 
O 

a^ 

> 

in 
o 

OO 

o\ 



X3 

o 



X 
S3 



We present a detailed account of the GW space-time 
method. The method increases the size of systems whose 
electronic structure can be studied with a computational im- 
plementation of Hedin's GW approximation. At the heart of 
the method is a representation of the Green's function G and 
the screened Coulomb interaction W in the real-space and 
imaginary-time domain, which allows a more efficient com- 
putation of the self-energy approximation E = iGW . For 
intermediate steps we freely change between representations 
in real and reciprocal space on the one hand, and imaginary 
time and imaginary energy on the other, using fast Fourier 
transforms. The power of the method is demonstrated using 
the example of Si with artificially increased unit cell sizes, 
keywords: electronic structure, quasiparticle energies, self- 
energy calculations, GW approximation 
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I. INTRODUCTION 

Computational electronic structure theory for real ma- 
terials depends on the use of simplifying approxima- 
tions for the many-electron problem. Two successful 
approaches have been the use of density functional the- 
ory and many-body perturbation theory. The density 
functional approach is overwhelmingly dominated by the 
local-density approximation (LDA) of Kohn and Sham 
[01, and extensions thereof, such as gradient corrections 
or self-interaction corrections. The many-electron prob- 
lem is mapped onto an effective non-interacting electron 
problem and solved for the ground-state density and en- 
ergy. The limitation of this approach lies in the fact that 
in principle it gives no access to the excitation spectrum 
of the system under study, even if approximations for the 
exchange and correlation potential are further refined. 
This limitation is felt particularly severely in semicon- 
ductor physics, where many of the phenomena of inter- 
est are centered on properties of the excited states. For 
this class of materials Hybertsen and Louie |0,^ showed 
that the GW approximation, first proposed by Hedin ^] 
in 1965, allows computation of band gaps in remarkably 
good agreement with experiment for a series of semicon- 
ducting and insulating materials. 

The GW approximation gives a comparatively sim- 
ple expression for the self-energy operator, which allows 
the one-particle Green's function of an interacting many- 
electron system to be described in terms of the Green's 
function of a hypothetical non-interacting system with 



an effective potential. The Green's function contains in- 
formation not only about the ground state density and 
energy but also about the quasiparticle spectrum. The 
GW approximation has still proved computationally very 
expensive and has mainly been used to determine the 
quasiparticle spectrum of bulk semiconductors and insu- 
lators although progress has also been made in 
the treatment of systems such as surfaces Q , clusters Q 
and simple polymers . The issue of self-consistency has 
only begun to be addressed fairly recently JlO|-p^. The 
situation is similar for total energy calculations, where 
thorough investigations so far have been made only for 
the homogeneous electron gas First GW calcula- 
tions of the charge density of Si and Gc have been per- 
formed recently urn. 

The method we describe in this paper substantially 
reduces the computational effort needed to study larger 
systems. The core idea was outlined in a Letter by Rojas, 
Godby, and Needs . Here we give a full account of the 
method and detailed aspects of our implementation. The 
improvements in efficiency over traditional implementa- 
tions of the GW approximation in a reciprocal-space for- 
malism result from choosing the representation most suit- 
able to the computational step being undertaken, either 
reciprocal space or real space on the one hand and imag- 
inary time or imaginary energy on the other, and switch- 
ing between representations easily with the help of fast 
Fourier transforms (FFTs). The choice of representing 
the time/energy dependence on the imaginary instead of 
on the real axis allows us to deal with smooth, decaying 
quantities which give faster convergence. To obtain the 
self-energy eventually on the real energy axis, we fit a 
model function to the computed self-energy on the imag- 
inary axis, and continue it analytically to the real axis. 

The paper is organised as follows: The second Section 
reviews briefly the equations needed for implementation 
of the GW approximation. Section III describes the es- 
sential concepts of the GW space-time method. In the 
fourth Section we outline detailed aspects of the imple- 
mentation such as mesh discretisation and exploitation 
of symmetry. In Section V we discuss the analytic con- 
tinuation procedure for the self-energy, and subsequent 
determination of the quasiparticle energies. In the sixth 
Section we present results for bulk silicon and silicon with 
artificially increased unit cell sizes, and discuss scaling 
behaviour for these examples. A summary concludes the 
paper. Atomic units are used throughout. 
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II. THE GW APPROXIMATION 

The central idea of Hedin's GW approximation is to 
approximate the self-energy operator E by 

C30 

S(r,r';c.)-^ J dc.' iy(r, r'; c.') G(r, r'; + c^')e'"'', 

— CO 

(2.1) 

where S is an infinitesimally small positive time, W is the 
screened Coulomb interaction, 



W{r,r';uj)^ J d^r"v{r ~ r") e-\r" ,r';Lj) (2.2) 

where v{r~r") is the Coulomb interaction l/|r — r"| and 
G is the one-particle Green's function. G itself depends 
on E through the Dyson equation and should arguably 
be determined self-consistently. In practice, however, in 
almost all calculations for real systems G has been ap- 
proximated by the non-interacting Green's function at 
the LDA level, i. e.. 



G^^-'(r,r';u;) 



nk 



v l/„k(r) 



(2.3) 



where 77 is a positive (negative) infinitesimal for occupied 
(unoccupied) one-particle states. The wavefunctions 'I'„k 
in this equation are eigenfunctions, with eigenvalues Snk, 
determined from a self-consistent LDA calculation for the 
system under consideration. 



For the inverse dielectric function in Eq. (2.2) one has 
to rely on suitable approximations. We use the random 
phase approximation (RPA), 



e^^^(r,r';a;) =5(r-r') 



(2.4) 



-J dv"v{v-v") P\v",v'-uo) 



with the irreducible polarization propagator P° at RPA 
level given by 

P"(r,r';a;) (2.5) 
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dw' G^^^(r, v';oj') G^^^(r, r'; w' - uj) 



Part of the efficiency of our method deriv es fr om the 
fact that the convolutions Eqs. (2.1) and (2.5) in the 
frequency domain become simple multiplications in the 
time domain (real or imaginary): s ee n ext section. For 
real times the Green's function Eq. (2_^) becomes 

G^^^(r,r';r) (2.6) 



i E *nk(r)*;k(i'') exp(-ie„kr), r < 0, 

nk 
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-i E ^'nk(r)*;k(i"')exp(-ie„kT), r > 0. 

nk 



For imaginary times our expression for G^^^ (see Eq. 
( |3.3|) later) corresponds to analytically continuing the 
T < form (the retarded Green's function) to the positive 
imaginary time axis, and the r > form (the advanced 
Green's function) to the negative imaginary axis. 

Once the self-energy operator is known we can employ 
first-order perturbation theory in — V^J^^) to com- 
pute quasiparticle corrections to the LDA eigenenergies 
(see Section VB). 



III. THE GW SPACE-TIME METHOD 
A. Mathematical formulation 

The traditional way to set up and solve the equations 
for the GW approximation has been to express and com- 
pute all quantities in the reciprocal-space and energy do- 
main, using a plasmon-pole model for the energy depen- 
dence of T4^ ||. To compute P° (Eq. ^) and S (Eq. 
(^.l| )) this involves sums scaling with the fourth power 
of the number of plane waves Nc used to represent the 
wavefunctions and quadratically with the number of en- 
ergy points used to represent the energy dependence, 
whereas the scaling in the real space and time domain 
is quadratic in Ng and linear in N^, since convolutions 
in the reciprocal-space and energy domain become sim- 
ple multiplications in the real-space and time domain. 
The expressions for e and W, on the other hand, are 
more efficiently calculated in a reciprocal-space and en- 
ergy representation. As fast Fourier transforms permit 
us to switch between representations with a numerical 
cost that scales like TV log TV, where iV is the number of 
points involved in the EFT, we can efficiently exploit the 
advantages offered by the one or other representation. 

The time- or energy-dependence of the quantities in- 
volved shows a structure that is not easily represented 
on an equidistant grid suitable for an EFT. However, 
we can rigorously analytically continue the quantities to 
the imaginary time or energy axes where the structure is 
much smoother and therefore amenable to a representa- 
tion on an equidistant grid. This point is illustrated by 
Fig. 1 of Ref. ||l6| which shows the energy dependence on 
the real and imaginary axis of the imaginary part of the 
self-energy Y^{k,uj) of jellium with a density parameter 
of Ts = 2.0. The rather ragged shape on the real energy 
axis contrasts with the smooth shape on the imaginary 
energy axis. We emphasise that, while the same amount 
of physical information is contained for E or similar func- 
tions on the real or imaginary time or energy axis, to 
obtain a well converged final answer for quantities that 
go through several Fourier transformations, the impor- 
tant information is more easily represented in imaginary 
time/energy, in much the same way that the choice of 
basis set can reduce the effort needed to satisfactorily 
represent wavefunctions. This is also shown in plots of 
the self-energy in Section VI, where it can be seen that 
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the smooth form of the self-energy on the imaginary axis 
stih aUows stable and accurate reproduction of the more 
complicated behaviour on the real axis. 

Using the imaginary time/energy representation allows 
us to explicitly take the time- or energy-dependence into 
account without having to rely on a plasmon-pole or 
other model for most of the calculation. Only after the 
full imaginary-energy dependence of the expectation val- 
ues of the self-energy operator has been established do we 
use a fitted model function (whose sophistication may be 
increased as necessary with negligible expense), which 
we then analytically continue to the real energy axis to 
compute the quasiparticle energies. 

The Fourier transforms between the complex axes work 
like their counterparts on the real axes, except that ad- 
ditional factors of ±i have to be included: 



F{it) = 



27r 



dujF{iLu) exp(zct;r). 



oo 

F{i(jj) — —i j dTF{iT) exp{—iujT) 



(3.1) 



(3.2) 



Mathematically they can be understood as Laplace trans- 
forms followed by analytic continuation to the imaginary 
axis. 

The computational steps which are successively under- 
taken in the GW space-time method are in detail: 

1. Construction of the Green's function in real space 
and imaginary time: pM 



G^^'4(r,rVr) 



(3.3) 



i ^nk(r)*;k(i"') exp(e„kT), r > 0, 



nk 
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-I E *„k(r)*;k(i-')exp(e„kT), r < 0, 

nk 



2. formation of the RPA irreducible polarizability in 
real space and imaginary time: 

F°(r,r';ir) = -iG^^'^(r, r'; iT)G^^'^(r', r; -ir), 

(3.4) 

3. Fourier transformation of to reciprocal space 
and imaginary energy and construction of the sym- 
metrised dielectric matrix in reciprocal space, 

e(k,G,G';zc^) =(5gg' 



4. inversion of the symmetrised dielectric matrix for 
each k point and each imaginary energy in recipro- 
cal space, 

5. calculation of the screened Coulomb interaction in 
reciprocal space: 



l¥(k,G,G';iw) 
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|k + G||k + G' 



X r\k,G,G';iw), (3.6) 

6. Fourier transformation of W to real space and 
imaginary time, 

7. computation of the self-energy operator: 

I](r, r'; ir) = iG(r, r'; ir)M^(r, r'; ir), (3.7) 

8. evaluation of the expectation values: 

(^'„k|S(ir)|*„k), (3.8) 

9. Fourier transformation of the expectation values to 
imaginary energy, 

10. fitting of a model function to the expectation val- 
ues of the self-energy, allowing analytic continua- 
tion onto the real energy axis, 

11. evaluation of the quasiparticle corrections to the 
LDA eigenvalues by first-order perturbation theory 
in (S-K^DA^. 

B. Discretisation of the equations 

A practical implementation on a computer requires the 
integrals described in the last section to be discretised 
and suitably truncated. Exploitation of symmetry can 
help to keep the computational cost for real materials 
down. These issues will be addressed in this section. 

The quantities we are dealing with in real space, such 
as G, P, VF, and S are non-local operators that decay 
as their two spatial arguments move apart. With the ex- 
ception of W ^ which will receive some extra attention, 
the non-locality is short-ranged, i.e, decaying faster than 
|r — r'l^^. This suggests that the distance the two vari- 
ables are allowed to move apart can be restricted suitably 
and the functions can be assumed zero beyond this range. 
We call this range the interaction cell (IC) . Furthermore, 
by the symmetry of the crystal, one of the arguments of 
any of these operators, which shall be symbolically de- 
noted here by P(r, r'), can be restricted to the irreducible 
wedge of one unit cell (lUC). The coarseness of the grid 
and the size of the interaction cell determine the pre- 
cision. A shape of interaction cell which is compatible 
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with the fast Fourier transforms and preserves symme- 
try at the same time is the Wigner-Seitz cell of a lattice 
whose defining vectors are multiples of the primitive vec- 
tors of the crystal lattice. We call this lattice the IC 
lattice (ICL) (see Figure |^) . We choose equal spacing for 
the lUC and IC grids. They could be offset from each 
other which would avoid the singularity of the Coulomb 
potential in real space but requires the handling of ad- 
ditional phase factors. We therefore normally choose the 
two grids to have no offset between them and deal with a 
single real-space grid (RSG), which is defined by vectors 
which are integer fractions of the primitive lattice vec- 
tors. The treatment of the sing ularit y of the C oulomb 
potential is discussed in sections IV B and IV C below. 

To summarise: the crystal is defined by the three prim- 
itive vectors ai, a2, a^, the ICL by the three vectors 



1, 



1,2,3, 



and the RSG by three vectors 

s,^SL,/N^, i = 1,2,3. 



(3.9) 



(3.10) 



In the example shown in Fig. |^ the Nf^ are 2 and the 
are 3. 

The grid vectors x in the IC are integer linear combi- 
nations of the Si and fulfil the Wigner-Seitz condition 



x-L < 

- 2 



(3.11) 



for any vector L that is an integer linear combination of 
the vectors 1. Only one of possibly several vectors for 
which the equality in Eq. (3.11) holds and which differ 
only by a lattice vector of the ICL must be contained in 
the IC, or, expressed differently, only half the surface of 
the Wigner-Seitz cell defined by Eq. (3.11) is part of the 
IC. The integers Nf^ in Eq. (3.9) determine the shape 



and size of the interaction cell in real space and the k- 
point grid in reciprocal space, as explained below. If 
there is no reason to believe that the non-locality of the 
operators has a very different range in one direction than 
in another, they should be chosen in such proportion to 
each other that the shape of the IC is as close to a sphere 
as possible. 

The first argument r of the functions of type F (see 
above) is again an integer linear combination of the vec- 
tors s and can be restricted to the irreducible wedge of 
one unit cell of the crystal lattice. This unit cell can in 
principle have any shape, but it helps to avoid having to 
deal with phase factors when applying symmetry opera- 
tions if this unit cell is also chosen to have Wigner-Seitz 
shape, i. e.. 



^ 1 9 
r • R < 



(3.12) 



where R is any crystal lattice vector and of any two r 
which are connected by a symmetry operation in the crys- 
tal's space group only one is kept. 



The choice of lattices in real space uniquely deter- 
mines conjugate lattices in reciprocal space and vice 
versa, if one follows the rule that the discrete approxi- 
mation to the Fourier transform must leave the functions 
unchanged if applied forwards and backwards in succes- 
sion. In reciprocal space, functions of type F (see above) 
are functions of the two variables G and K, where G 
is a reciprocal lattice vector and K can be written as 
K = k + G' , where G' is again a reciprocal lattice vec- 
tor and k a reciprocal vector in the first Brillouin Zone 
(BZ). K is restricted to the Wigner-Seitz cell of the re- 
ciprocal lattice of the RSG, which we shall call simply 
the reciprocal- space cell (RC) and G is restricted to its 
irreducible wedge (IRC). The spacing of the k is deter- 
mined by the primitive vectors of the reciprocal ICL and 
we shall call the grid that is defined that way the K grid 
(KG) ||. 

The transformation from real-space to reciprocal-space 
representation of F is given by 



F(G,K) = 

lUC star(r) IC(t) 



(3.13) 



r p r' 



p denotes a point group element (including possibly a 
non-primitive translation in non-symmorphic symmetry 
groups) and the sum over p runs over all symmetry op- 
erations that generate the full star of r, i. e., the set of 
all RSG points in the Wigner-Seitz cell that are obtained 
by applying symmetry operations of the crystal to this 
particular RSG point r in the lUC. Njc is the number 
of RSG points contained in the IC. The reverse transfor- 
mation is given by 

F(r,r')= (3.14) 

^ IRC star{G) RC 

y yi^(G,K)e'''-°e-^''("-'-"-)-^. 



Nuc 



G 



K 



The sum over p runs here over all symmetry operations 
that generate the full star of G and Njjc is the number 
of RSG points contained in one unit cell. 

Because of the symmetry of the problem and the fact 
that K = k -I- G', F can in reciprocal space be alter- 
natively written as a square matrix i^k(G,G') with k 
restricted to the irreducible wedge of the Brillouin Zone 
and G and G' given everywhere in the RC. k is here 
written as an index to emphasise the matrix nature of F 
in this representation. 

It should be briefly mentioned that the 'natural' cell 
form for a three-dimensional FFT is a parallelepiped and 
not a Wigner-Seitz cell such as we are dealing with. How- 
ever, there is a unique mapping between the two shapes 
using translations by vectors of the lattice defining the 
Wigner-Seitz cell in question, because F is invariant un- 
der such translations by construction. 
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Apart from the real-space and reciprocal-space repre- 
sentations we will occasionally use a mixed-space repre- 
sentation. This is defined by 



i^k(r,r') 



IC 

E 

R 



i^R(r,r')exp(-ik-R), 



(3.15) 



with the reverse transformation 

FR(r, r') = -r^ V Fk(r, r') exp (zk • R). (3.16) 
A/r ^ 

When we use the mixed-space representation, r and r' 
are always understood to be confined to the UC (or its 
irreducible wedge, see below). The notation FR(r, r') is 
another way of writing Fir -f R, r') and has been chosen 
to emphasise that any point on the RSG in the IC can 
be written as the sum of a crystal la ttice vector R and 
a vector in the UC. The sum in Eq. ( 3.15| ) runs over all 
crystal lattice vectors in the IC. The k are vectors in the 
first BZ of the crystal and are the conjugate vectors of 
the R. Symmetry can be exploited to reduce the number 
of points needed to represent F in mixed space. For 
example, r can be restricted to the irreducible wedge of 
the UC while r' is given for every point of the RSG in 
the UC and k on every point of the KG within the first 
BZ, or k can be restricted to the irreducible wedge of the 
first BZ, with r then given on every RSG point in the 
UC. 

A mixed-space (MS) formalism was used by Blase et 
al. for calculating the polarizability and dielectric 
function of periodic systems. A mixed-space representa- 
tion was also used and described by Godby, Sham and 
Schliiter |21| for GW calculations, albeit without refer- 
ring to it by that name. Blase et al. |Q discuss the 
computational efficiency of the MS scheme in compari- 
son with a direct real-space (RS) approach. In Fig. 2 
of Ref. they show how the number of (r, r') pairs 
which have to be computed in order to set up the polar- 
izability P^(r, r') for one k in the MS method (using a 
4x4x4 Monkhorst-Pack k grid) depends on the range 
of interaction (non- locality range) Rmax- They compare 
this to the number of (r, r') pairs obtained by confining 
jr — r'l to a sphere of radius Rmax which is assumed to 
be the corresponding number of pairs to be computed in 
a RS method. This is somewhat misleading since - as is 
discussed above - the size of the IC in real space is deter- 
mined by the size of the k grid of the conjugate mixed- or 
reciprocal-space quantity. Hence it does not make sense 
to increase the interaction sphere beyond the boundaries 
of this IC which correspond to a non-locality range of 
Rmax—^^-^ a.u. for a 4 X 4 X 4 k grid. Furthermore, the 
real-space quantity P°(r, r') contains information on all 
k points. So for a proper comparison between the two 
methods the number of pairs given in Ref. for the 
MS approach has to be multiplied by the number of spe- 
cial k points which is 10 for a 4 x 4 x 4 Monkhorst-Pack 
grid. The statement made in Ref. that the calcula- 
tion of P*'(r,r',cj) requires a double BZ summation in 



the real-space scheme - as opposed to a single BZ sum- 
mation in the MS scheme - does not hold any more if 
the Green's function is set up in mixed space and then 
Fourier transformed to real space, as described in section 
IV Al below. 



So far we have established the nature of the grids in 
real and reciprocal space that are suitable for use with 
the discrete Fourier transform. These grids fill a Wigner- 
Seitz-cell shaped volume which in turn can be uniquely 
mapped onto a parallelepiped, the shape that is required 
by the discrete Fourier transforms. However, the starting 
point of the GW calculation is the output of a standard 
LDA calculation using plane waves and pseudopotentials 
which determines the Fourier coefficients of the wave- 
functions ^'„k(G) for all reciprocal lattice vectors that 
lie within a sphere defined by (k -f G)^ less than some 
cutoff. The RC must therefore be big enough to comprise 
all of these 'shifted spheres'. The volume of the small- 
est possible such cell will be typically 2 to 4 times larger 
than the volume of the individual spheres. Furthermore, 
to strictly prevent aliasing in the steps where we replace 
a convolution in one representation by a multiplication 
in the conjugate representation, we would have to choose 
the grid for the FFTs twice as large as this minimal cell 
in every dimension, increasing the volume by a factor 
of eight. This means that the FFT grid would have to 
comprise between 16 and 32 times more points than the 
initial LDA calculation had reciprocal lattice vectors. We 
have found that in practice aliasing effects are very small 
once we choose enough plane waves for good overall con- 
vergence, and in the limit of an infinite number of plane 
waves any aliasing effects vanish strictly along with any 
other error caused by truncation of the set of recipro- 
cal lattice vectors. Therefore it would not be justified 
to accept the huge overhead imposed by the doubling of 
the grid size in all dimensions, and we usually restrict 
our FFT grid to the minimum Wigner-Seitz cell as de- 
fined above, leaving us with an FFT grid with between 2 
and 4 times as many points as the LDA calculation used 
reciprocal lattice vectors. 

The computational effort can be further reduced by 
physical considerations. The FFT grid fills a Wigner- 
Seitz cell shaped volume in either space. However, if 
we assume in real space that the range of non-locality 
is uniform in all directions we can neglect all the points 
in the interaction cell outside the largest possible sphere 
inscribed into it. Similarly we can usually assume that 
the Fourier coefficients in reciprocal space fall to zero 
after an equal length in all directions and thus we can 
cut back in reciprocal space to the largest possible sphere 
inscribed into the RC lattice. It has to be kept in mind 
that for symmetry reasons it is the vectors k -t- G which 
must fit into the sphere and not merely the vectors G. 
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IV. NUMERICAL ASPECTS AND SCALING 
WITH SYSTEM SIZE 

As the GW space-time method is primarily designed 
to enable larger systems to be studied within the frame- 
work of the GW approximation, in this section we study 
the scaling of the computational cost with system size. 
To gauge the cost of a calculation, we first look at the 
real-space, imaginary-time representation. What mat- 
ters here is the number of points in the irreducible wedge 
of the unit cell, the number of points in the interaction 
cell and the number of points on the imaginary time axis. 
We will discuss the time dependence of the operators in 
the next section and concentrate here on the spatial di- 
mensions. Since in our setup the grid spacing in the 
unit cell and in the interaction cell are chosen equal, the 
three factors determining the cost of the calculation are 
the size of the lUC, which is a system property, and the 
size of the interaction cell and the grid spacing, which 
are convergence parameters. 



A. Green's function and polarizability 

The setting up of the Green's function is the starting 
point for the GW calculation. We take the Fourier coeffi- 
cients of the lattice-periodic part of the Bloch functions, 
Mnk(G), and the eigenvalues e„k from a previous stan- 
dard LDA calculation. The eigenvalues are expressed on 
an energy scale with zero at the Fermi energy, which for 
semiconductors and insulators is taken to be in the mid- 
dle of the band gap. The wavefunctions u„k are trans- 
formed via an FFT to real space and the Green's function 
initially computed in mixed space: 



R 



G£^■'(r,r';^r) 



(4.1) 



= < 



I J2 *nk(r)*;k(i'') exp(e„kT), r > 0, 



E ^nkWl'^kli'Oexple^.kT), r<0. 



The number of unoccupied bands included in the 
Green's function for r < is a convergence parameter. 
Because of the rapid decay of the exponentials for higher 
energies it suffices to set a cutoff energy for bands in- 
cluded in the sum that is considerably smaller than the 
largest LDA eigenvalue, but this cutoff must be several 
times larger than the highest quasi-particle energy to be 
computed. We denote the total number of bands, occu- 
pied and unoccupied, included in the sum as Nbands- 

The Green's function can then be transformed from 
mixed space to real space using an FFT: 



ir) exp(ik • R), 
(4.2) 



where the sum goes over the k grid in the Brillouin zone 
corresponding to the interaction cell (see III(B) ear lier). 

To set up the function in mixed space (Eq. (4.1)) we 
need Njuc x Njjc x Nbands multiplications in the spatial 
dimensions. The transformation to real space brings in 
another factor of Nn, the number of unit cells contained 
in the interaction cell, (disregarding the additional fac- 
tor of logiVR. from the FFT), so that the overall scal- 
ing is like Niuc x Nic x Nbands- Since the number of 
bands is determined by a cutoff energy that is constant 
for a given material, Nbands will grow linearly with sys- 
tem size. The computational cost at this stage therefore 
scales quadratically with system size assuming a fixed 
size of interaction cell (i. e., fixed range of non- locality). 
Additional care is necessary once the unit cell outgrows 
the range of non-locality. In this case there is only a single 
k point necessary and the transformation between mixed 
and real space is redundant. To prevent a crossover to 
cubic scaling the interaction cell must be kept constant 
at its converged size and will then not be a multiple of 
the unit cell size, but smaller than it. 

The next stage, where the irreducible pol arizability is 
formed in real space according to Eq. ( ^.4[ ), scales lin- 
early. 



B. Dielectric matrix 

The formula for the dielectric function in real space 
and on the imaginary energy axis reads 

IC(r') 

e{r,r';iLu) = S{r-r')- v{r - r")P°{r" ,r';iu;). 

r" 

(4.3) 

Here r' is restricted to the lUC, r to within an IC around 
r', and r" in the sum on the right hand side runs over an 
interaction cell around r'. The right hand side contains 
a convolution that can be more efficiently handled in re- 
ciprocal space where it becomes a simple multiplication. 
As only one of the two spatial arguments of the dielectric 
function is involved, we can write 



e(K,r') = t.(K)P"(K,r'), 



(4.4) 



where K is the conjugate variable of x = r' — r in recipro- 
cal space. (The imaginary-energy argument is suppressed 
here and until the end of this subsection.) Again we can 
see that for a given interaction cell size the computational 
effort to set up the dielectric function scales linearly with 
the number of points in the lUC. 

The long-range behaviour of the diele ctri c function re- 
quires some special consideration. Eq. (O) contains the 
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convolution of the long-ranged Coulomb potential with 
the short-ranged non-locality of P°. In fact, for any r' 



J drP"(r,r') 



0, 



whereas 



drv{r) 



(4.5) 



(4.6) 



The analytic convolution of the two functions over all 
space yields a function whose integral over all space is 
finite, meaning that the Fourier coefficient for K = 
has a finite, non-zero value. To preserve this behaviour, 
which is important if one is actually interested in comput- 
ing the value of the dielectric constant but also for fast 
convergence with respect to interaction cell size of the 
quasiparticle energi es, o ne has to use a modified Coulomb 
interaction in Eg. (|4.3|) . By construction, the property 
expressed by Eq. (|4.5D is within our discrete and strictly 
finite range approximation compressed into the interac- 
tion cell: 



/C(r': 

E 



P°(r,r') = 0. 



(4.7) 



To obtain the right dielectric constan t, w e have to com- 
press the property expressed by Eq. (4.6) into the finite 
interaction cell as well. A natural way to do this is to 
use the reciprocal-space definition of v: 



47r 



(4,8) 



Transforming to real space, through a discrete Fourier 
transform restricted to one interaction cell, yields an ex- 
pression that contains correction terms to the simple 1 /r 
form that will vanish in the limit of infinitesimal grid 
spacing and infinite IC s ize. Since, in practice, we evalu- 
ate the convolution Eq. (4.3) as a multiplication in recip- 
rocal space, there is no need to explicitly transform Eq. 
(18) to real space. 

In practice, to deal with the divergence of the Coulomb 
potential at zero wavevectors, we follow the procedure 
described in the literature [^2yi8[ , i. e., we employ k ■ p 
perturbation theory for calculating the head (k = G = 
G' = 0) and wings (k = G = 0; G' ^ 0) of the sym- 
metrised dielectric matrix Eq. (3.5). The matrix elements 



of type (4'nq|e~ |^n,q-k) appearing in the expression 
for the RPA polarizability in reciprocal-space represen- 
tation (cf. Eq. (18) of Ref. [^) are evaluated by expand- 
ing the wavevector dependence of the wavefunctions and 
the exponential for small k. This yields the lowest-order 
(in k) terms for head and wings of the polarizability. 
The head of the polarizability goes to zero like k'^ thus 
cancelling the divergence of the Coulomb poten- 

tial at k = G = G' = 0, yielding a finite value for the 
head of uP" whereas the lowest-order term for the wings 



Pqq, (0) is proportional to k, cancelling the 1/k diver- 
gence of woG'(O) and yielding a finite (t;P'')oG' (0). Head, 
wings and body of the inverse dielectric matrix are then 
obtained in terms of head, wings and body of the dielec- 
tric matrix by block-wise inversion as described by Pick, 
Cohen and Martin [||. 

The inversion of the dielectric matrix can be performed 
either in mixed space or in reciprocal space. While this 
is strictly equivalent if the whole FFT grid is used in ei- 
ther space, it is important to look at ways to reduce the 
grid size without compromising the exactness of the re- 
sult, as the matrix inversion is in fact the one operation 
in the computational sequence which scales worst with 
system size (unless sparseness is exploited, see below) 
and will therefore be the bottleneck for large systems. In 
mixed space we have to invert square matrices of dimen- 
sion Nuc for every k point in the irreducible wedge of 
the Brillouin zone (IBZ). 

First we look at the case where the IC is in no dimen- 
sion smaller than the UC. We are then dealing with the 
inversion of fully occupied square matrices, which scales 
cubicly with the dimension of the matrix, so that the in- 
versions together scale like Njbz >^ -^uc- This number 
is essentially the same as Nujc x -^/c x -^c/C, small dif- 
ferences arising possibly because of the coarseness of the 
grid. This shows that the scaling is quadratic with sys- 
tem size for constant IC size and constant Njjc / Nijjc ■ 
If the increase of the unit cell size means also a lowering 
of symmetry, an additional factor corresponding to the 
symmetry reduction comes in. Transforming to recipro- 
cal space and cutting back to a set of reciprocal lattice 
vectors G within a sphere as described earlier reduces 
the dimension of the matrices by a factor of 2 to 4, which 
speeds up the inversions by a factor of 8 to 64, because 
of the cubic scaling. This is the procedure we choose 
routinely in our calculations. 

If the unit cell size exceeds the interaction cell size 
in one or more dimensions the matrix in mixed space be- 
comes sparse, as all those elements for which r' is outside 
the IC become zero. This sparseness can be exploited to 
the effect that the matrix inversion scales only as A^^c", 
so that the overall scaling with system size remains the 
same. However, in the scheme we are employing at the 
moment this sparseness is not yet explicitly exploited so 
that we are in fact dealing with a situation where the IC 
grows with the UC, leading for the matrix inversion to an 
Nfj(j scaling once the unit cell outgrows the interaction 
cell. 



C. Screened Coulomb interaction 

The screened Coulomb interaction is defined in real 
space as 

W{Y,Y']iiu) = I dr"e-^{r,r";iuj)v{r" -r') (4.9) 
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The convolution on the right hand side is again best dealt 
with as a multiplication in reciprocal space 

Wgg' (k; iuj) = ecG' (k; iuj)vGG' (k). (4.10) 

However, in a semiconductor and in metals for frequen- 
cies other than zero, W is truly long ranged, and there- 
fore diverges for zero wavevectors. In order to avoid prob- 
lems (in the G/r and G'/r' FFT) associated with the re- 
sulting long-range tail and (in the It jiu) FFT) with the 
asymptotic frequency dependence we separate W into a 
long-range part which is immediately set up in real space 
and a short-range part Ws which is first computed in 
reciprocal space. From the long-range behavior of the di- 
electric function we find the isotropic part of the screened 
interaction in the long-range limit 

lim VK(r,r';it^) = /(it^)w(r-r') (4.11) 

|r— r' I — >oc 

SO that we can define a short-ranged part Ws of W 



W,(r,r';zw) = 




r,r";zt^)-/(itj)(5(r-r")) 



xw(r"-r'). (4.12) 

Ws has well-defined Fourier coefficients for all recipro- 
cal vectors. It can thus be computed in reciprocal space. 
The long-range part is set up in real space with /(iw) 
determined from the small- wavevector behavior of the in- 
verse dielectric matrix in reciprocal space, i. e., the head 
of the inverse dielectric matrix e,^^: 

(4.13) 

Even though the interaction is truly long-ranged we 
only need to set it up inside the IC because we aim to 
calculate the self-energy whose range is determined by 
that of the (short-ranged) Green's function Gq. Hence, 
at this point, it is sufficient to take the contribution of 
the long-range part of W into account properly within the 
IC. The Coulomb potential at r = r' is approximated by 
the average over a sphere around r = r' with the same 
volume as one real-space grid cell. 

D. Self-energy 

The self-energy operator is set up in the real-space and 
imaginary-time domain: 

S(r, r'; ir) = iG^^^(r, r'; iT)W{r, r'; ir - i6). (4.14) 

Within the discrete-grid and finite-range approximation 
a complete description is given if r is restricted to the lUC 
and r' to an IC around r. The scaling is again linear with 
Njuc for fixed Njc- 



E. Expectation values of the self-energy 

To compute the expectation values of E between wave- 
functions at the special k points, we transform E to 
mixed space and form the matrix elements 

uc 

(vI/„fclE(zu;)|vI'„fc} = ^vl/;^(r)Ek(r,r';zc^)vI/„k(r') 

r.r' 

(4.15) 

This operation again scales quadratically with the num- 
ber of points in the UC for any given set of quantum 
numbers nk. For a general point q in the first Brillouin 
zone that is not on our discrete grid, we generate inter- 
polated values through the formula 

IC 

Sq = ^ERexp(-zq-R). (4.16) 

R 



V. QUASIPARTICLE ENERGIES 
A. Analytic continuation 

In order to calculate the quasiparticle corrections to 
the LDA eigenvalues we require the self-energy operator 
as a function of real energy, and thus a key requirement 
for our imaginary time method is the ability to obtain 
the expectation values of E on the real energy axis accu- 
rately from the imaginary energy behaviour. From com- 
plex analysis we know that if two functions are equal 
over any arc in the complex plane then they are equal 
everywhere in their common region of analyticity. We 
know from the structure of G and W that T,{z) {z denot- 
ing complex energy) has poles in the second and fourth 
quadrant of the complex plane. If we know the ana- 
lytic form of the expectation values of the self-energy on 
the imaginary energy axis we can analytically continue it 
from the negative imaginary energy axis to the negative 
real energy axis, and from the positive imaginary to the 
positive real axis, without crossing any branch cuts. 

To obtain such an analytic form, we fit a multipole 
model function for each pair of quantum numbers nk to 
the values (*„k|E(icj)|^'„k) = {^(i^))- 

71 I 

(vI'„k|E(^)|*„k) ^a°k + E^^' (5-1) 

where a^^. and 6*^^. are complex fit parameters and z is 
the complex energy. In Fig. |^ the resulting fitted form of 
the correlation (energy-dependent) part of the self-energy 
is compared with the calculated matrix elements for 
bands at the F and X points (the fits are plotted with the 
same line styles as the respective calculated matrix ele- 
ments, they are on this scale indistinguishable from the 
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latter). A simple two-pole model (n = 2) performs very 
successfully for Si, with the fitted function reproducing 
the actual values with an rms error of less than 0.2%. 
Including several further poles in the fitted form proves 
to be stable but unnecessary, although more poles are ex- 
pected to be required for systems with multiple natural 
energy scales. For the analytic continuation to be valid, 
the fitted pole positions 6^^. should lie in the upper half 
plane when fitting the negative imaginary axis and vice 
versa, a condition which in practice is obeyed by the opti- 
mal parameters. The parameters a*^^ should in principle 
be zero, as limi^^oo i^i'-^)) = 0- Allowing a small finite 
value for a°i^ has proved helpful in the fits, though, and 
as we are interested mainly in energies close to the Fermi 
energy and not the long-range limit this is perfectly le- 
gitimate. Since = iw))* it is sufficient to fit 
only one half-axis. 

Convergence of the quasiparticle energies with respect 
to the parameters of the time-energy transform grid is 
discussed in the next section, but it is also important to 
consider as a separate question the stability of the ana- 
lytic continuation procedure. In particular if we are to 
achieve smooth convergence it is important that changes 
in the calculated self-energy on the real axis correspond 
to genuine changes in the calculated (E(ia;)), and are 
not simply resulting from instabilities in the fitting pro- 
cedure. In Figures I and I we show the convergence 
behaviour of a matrix element (S(ia;)) at F with respect 
to number and spacing of the energy points respectively, 
and the corresponding form of the analytically contin- 
ued element. It can be seen that the fitting approach 
is indeed stable, with the convergence of (S(a;)) being 
dictated directly by changes in the calculated (E(iti;)). 
Thus the task of converging the quasiparticle energies is 
equivalent to converging the self-energy on the imaginary 
axis, with very little comparative loss of accuracy in the 
analytic continuation procedure itself. 

The parameter tOmax (describing the energy range used 
in the calculation) mainly determines the convergence 
of the large-energy region of /m(Ec(*w)). The energy 
grid spacing Aw predominantly affects the convergence 
of (Ec(ici;)) in the region close to uj= 0. The behavior of 
(Sc(a;)) (on the whole real energy axis) is more sensitive 
to the shape of (Ediw)) in this region of the imaginary 
axis than to the large-imaginary-energy tail of the latter 
(see also next paragraph). That is why the convergence 
with respect to Aw shown in Fig. || appears to be better 
on the imaginary axis than on the real axis. 

Our investigations of the convergence of the calcu- 
lated matrix elements of the correlation part of the self- 
energy and the resulting quasiparticle energies with re- 
spect to ujmax have shown that good convergence can 
be achieved by keeping the energy range for the fitting 
procedure fixed rather than fitting the whole range of en- 
ergies {—LUmax, +^max) whcu increasing the latter. The 
energy range for fitting the matrix elements {Y,c{iuj)) was 
fixed to (—5 Har, -1-5 Har) in the present calculation (the 
results are not sensitive to the exact value). The rea- 



son for restricting the fit range for (Ec(ia;)) is that most 
information is contained in the form of this function at 
imaginary energies reasonably close to iut = 0. Thus, fit- 
ting a large range, beyond a certain point, will mean that 
this part of the shape will be less accurately described due 
to the loss of weight in this region. This effect prevents 
(Ec(w)) from converging properly unless the fit range is 
kept fixed (at any reasonably sensible value), as is illus- 
trated by Fig. ^. 



B. Quasiparticle corrections 

The quasiparticle energies are formally solutions of the 
equations 



-^V2 + K.t(r) + V^H(r))vI/fJr) 



+ |dr'E(r,r';<f,)*rk(r') = <WW- (5-2) 

Because E is in general a non-Hermitian operator, the 
eigenenergies e^^ are complex, their real part being inter- 
preted as quasiparticle energy and their imaginary part 
as lifetime. It has been observed ||3,|2l|l that the wave- 



functions obeying Eq. (5.2) for typical semiconductors 
and insulators have an almost complete overlap with the 
LDA eigenf unctions which solve an equation in which 
the non-local self-energy operator is replaced by a local 
exchange-correlation potential. This makes it possible 
to find the quasiparticle energies by computing correc- 
tions to the LDA eigenvalues in first-order perturbation 
theory in most cases. (For some systems, however, the 
quasiparticle wavefunctions are expected to be qualita- 
tively different from the LDA wavefunctions, such as at a 
surface. The space-time method allows the full quasipar- 
ticle wavefunctions and energies to be calculated where 
required |2^.) Usually the LDA Hamiltonian with its 
eigenf unctions and eigenvalues is taken as the starting 
approximation, with E(w) — V^c as a perturbation. Fol- 
lowing an idea used by Hedin Q for the electron gas, we 
shift the LDA eigenergies by a constant that aligns the 
Fermi energies of the quasiparticles before and after ap- 
plying the GW correction. This is intended to simulate 
to some extent the effect of self-consistency in G and has 
been shown in a model system [ p6| to be instrumental in 
keeping charge conservation violations negligible. 

Calculation of the full energy dependence of the self- 
energy (via the analytic continuation to real energies) 
allows us to solve the equation 



,9P _ ,LDA 



(^-^.kl E(efk) - V^c - es l^-nk) (5.3) 



for the quasiparticle energies self-consistently. This ap- 
proach was employed for the calculation of the quasi- 
particle energies in the present paper. Alternatively, a 
Taylor expansion of E(w) at w = ^nW^ '^^^ ^le used: 
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,9P _ ^LDA 



I*nk) , 



duj 



(^'„k| l^'^.k) 



(5.4) 



(5.5) 



Because of the analytic fit of the expectation values of the 
self-energy the derivative with respect to the energy argu- 
ment is readily found in closed form. For eg a closed form 
but rather involved expression can also be given. Com- 
parison of the self-consistent solution o f E g. (5.3) and 



the result of the Taylor expansion Eq. (5.4) shows that 
the latter is usually a very good approximation, yielding 
quasiparticle energi es di ffering by only up to 5 meV from 
the solution of Eq. (5.3) for bulk Si. 



VI. 



RESULTS FOR BULK SILICON AND 
SILICON SUPERCELLS 



In order to test the performance and the scaling prop- 
erties of the space-time GW method we have performed 
calculations of self-energy and quasiparticle energies for 
bulk silicon and silicon with artificially increased unit cell 
sizes. The latter were slab-like supercells containing four 
and eight atoms, respectively, with tetragonal symmetry, 
i. e., the unit cells were artificially enlarged in one direc- 
tion. Supercells of this geometry can be used to model 
semiconductor superlattices and surfaces. 



A. Convergence parameters 

Before discussing the performance and scaling issues 
we briefly summarize the results of convergence tests for 
bulk Si. In Table | the parameters needed to converge 
quasiparticle energy differences to an accuracy of 50 meV 
and 20 meV, respectively, with the present method are 
gathered. This accuracy is meant to be with respect to 
each parameter individually. We estimate the overall ac- 
curacy of our results to be better than 0.1 eV. The plane- 
wave (PW) cutoff in Table | is the cutoff of i(k G)^ 
in the LDA calculation providing the wavefunctions ^'„k- 
All PW components of the LDA wave functions are used 
in the GW calc ulatio n. The real-space grids in the unit 
cell (see Section III B ) resulting from PW cutoffs of 13.5 
Ry and 16 Ry comprise 9x9x9 and 12 x 12 x 12 points, 
respectively. The band cutoff determines the nu mbe r of 
bands included in the band summation in Eq. ( [3.5| ) for 
the calculation of the Green's function. Band cutoffs of 
8 and 10 Ry correspond to taking 109 and 145 bands, re- 
spectively, into account. For sampling the BZ we employ 
a regular non-offset grid of k points. The LDA wave- 
functions are calculated at the k points in the irreducible 
wedge of the BZ. As described in Section IIIB the dimen- 
sions of the k grid in the BZ are equal to the dimensions 



(in multiples of primitive lattice vectors) Nj^ of the IC, 
which must be large enough to comprise the range of non- 
locality of the Green's function and the self-energy in Si. 
The smaller k grid given in Table | can roughly accomo- 
date a non-locality range of 14.5 a.u. and the larger of 
18 a.u. {—T„-iax, +Tmax) and At stand for the range and 
spacing of the equidistant grid for sampling the functions 
-F(r, r'jir) on the imaginary time axis. The parameters 
^max and Aw of the corresponding grid in the imaginary 
energy domain are related to the time-grid parameters 

by Aw = 2TT/{2Nr - 1)At with Tmax = NrAr. 



B. Quasiparticle energies for bulk silicon 

The quasiparticle energies for bulk Si obtained from 
a well converged calculation with the present method 
are given m Table [| The LDA calculation was 

done with the PW pseudopotential method employ- 
ing a pseudopotential constructed according to the pre- 
scription of Hamann and using the Ceperley- Alder 
exchange-correlation potential | 29 [| in the parametriza- 
tion of Perdew and Zunger ||30|. As can be seen from 
Table || the calculated GW quasiparticle energies agree 
well with quasiparticle excitation energies derived from 
photoemission, inverse photoemission and optical exper- 
iments pl|-|37t , except for the position of the bottom of 
the valence band which appears to be too high in our cal- 
culation. The agreement of the quasiparticle energy dif- 
ferences obtained with the present method {GW 1 in Ta- 
ble H) with those of a recent GW calculation by Fleszar 
and Hanke j38| is very good. Like us, these authors did 
not make use of a plasmon pole model to describe the 
energy dependence of the dynamically screened Coulomb 
interaction but directly computed its full energy depen- 
dence. 

We investigated how the quasiparticle energies change 
when the Green's function is updated by replacing the 
LDA eigenvalues by the quasiparticle energies and the 
self-energy recalculated with the updated Green's func- 
tion as input, as was done by Hybertsen and Louie 
The results {GW 2) are also included in Table ||. Updat- 
ing the Green's function increases the gap by about 0.1 
eV. Now the conduction band energies agree within bet- 
ter than 0.1 eV with those given in Ref. |^ whereas the 
valence band energies remain higher by between 0.1 eV 
and 0.5 eV (relative to the top of the valence band), the 
difference increasing with the distance of the respective 
state from the Fermi level. 

Indeed, the results of several earlier GW calculations 
Pl,^ are somewhat different from those given in Table || 
(and Ref. [Q). This can be attributed to the following 
reasons: 

(i) a possible lack of convergence with respect to the 
number of conduction bands taken into account in 
the earlier calculations, causing the topmost occu- 
pied state (which is particularly sensitive to this 
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parameter) to be about 0.2 eV too high, as dis- 
cussed in Ref. Q ; and 

(ii) the use of plasmon pole models for the energy de- 
pendence of the dynamically screened Coulomb in- 
teraction in the earlier calculations. 



C. Silicon supercells 



In Table IV we compare the quasiparticle energies 
for slab-like Si4 and Sis supercells to the bulk Si re- 
sults. Shown are the quasiparticle energies for the high- 
symmetry k points of Si (which are not necessarily sym- 
metry points for the supercell geometry) and for those 
k points which are mapped onto these symmetry points 
when the unit cell size is increased. The calculations were 
done with a PW cutoff of 13.5 Ry (leading to a real-space 
grid spacing of Ar = 0.8 a.u.), a band cutoff of 10 Ry, At 
= 0.3 a.u. and Tmax = 20 a.u. The k grids have been cho- 
sen to make the IC of equal size in all three dimensions 
for all the calculations, for the reason discussed in Sec- 
tion III B. The grid sizes resulting from these parameters 

and Nic are 



are summarized in Table III. N 



uuc 



uc, 



the numbers of real-space grid points in the irreducible 
part of the unit cell, the unit cell and the interaction 
cell, respectively (see Section IV above). Nk stands for 
the number of k points in the BZ and Nt/c^ is the number 
of positive time/energy points. The supercell results in 
Table agree with the bulk Si results, as expected. The 
small differences for the highest conduction band given 
m Table are ascribed to the lower symmetry of the 
supercells leading to different G vector sets. 



D. Scaling 



As outlined in Section IV the computational effort in 



the present method should scale quadratically with the 
number of atoms in the unit cell. Fig. ^ shows the CPU 
times for a full GW calculation for Si„ {n — 2, 4, 8) on 
a Digital Alpha 500/500 workstation. The parameters of 
these calculations are those given in Table III and the 
discussion thereof. Si2 (bulk Si) has a higher symme- 
try than the tetragonal supercells Si4 and Sis. As the 
symmetry is exploited in most parts of the calculation 
this saves approximately a factor of two in CPU time 
(cf. Niuc values in Table III). Taking this into account 
the overall scaling is indeed quadratic, as expected. The 
most time consuming parts of the calculation are: 

(i) setting up the Green's function according to Eq. 

(PD 



(ii) calculating the dynamically screened Coulomb in- 
teraction_(including the inversion of the dielectric 
matrix pq] , the transformation to and from recip- 
rocal space, and reading of P'^ from and writing of 
W to disk) 



(iii) computation of the matrix elements of the en- 
ergy dependent self-energy for a given number of 
k points and bands. 

Optionally, the Green's function can be recomputed when 
the self-energy is calculated, rather than storing it on disk 
when it is first set up and reading it in for the calculation 
of the self-energy. This procedure reduces the required 
disk space by a factor of 2/3. It has been used for the 
Si supercell calculations described here. As can be seen 
from the breakdown of the total CPU time shown in Fig. 
^ all these parts of the calculation scale roughly like N"^. 
More precisely, the scaling is like NnjcNbands for (i), 
NiucNuc for (h) and N^q for (iii) as was discussed in 



Section IV 



Comparing the performance of our method with con- 
ventional techniques we note that although a number of 
quasiparticle calculations within the GW approximation 
have been reported for systems with up to 60 atoms, 
e. g. for surfaces defects and fuUerenes Q, 

several simplifying approximations have been employed 
in these calculations in order to reduce the computa- 
tional effort. These authors focus on the calculation 
of quasiparticle energies and therefore only a number 
of self-energy matrix elements were computed. This is 
in contrast to the present method which provides the 
full self-energy (thus giving acccess to quantities other 
than the quasiparticle energies, see conclusions section 
below). In all the calculations mentioned above plasmon- 
pole models have been employed to approximate the en- 
ergy dependence of the dynamically screened Coulomb 
interaction. In some cases pT| , p4[ a model static dielec- 
tric matrix has been used as well. To aid comparison 
we note that although the plane-wave cutoff of the un- 
derlying LDA pseudopotential calculation can be (and 
has been in many GW calculations reported in the lit- 
erature) considerably reduced in the GW calculation be- 
cause the quasiparticle corrections usually converge much 
faster with this parameter than the LDA eigenvalues, this 
has not been done in the test calculations reported in the 
present paper. The energy dependence of the screened 
interaction is fully taken into account in the method pro- 
posed in this paper. Hence the dielectric matrix is com- 
puted and inverted for about 60 to 100 imaginary en- 
ergies which is computationally much more demanding 
than using a plasmon-pole model. We estimate that the 
crossover to the advantegeous scaling behavior of CPU 
time in our method occurs (a) already for bulk Si if we 
compare with a calculation of the full self-energy with a 
conventional reciprocal-space approach and (b) for sys- 
tems in the range studied in this paper if we compare with 
a reciprocal-space method where a plasmon-pole model 
is used and which is restricted to computing a moderate 
number of self-energy matrix elements only. Finally, we 
mention that work, particularly regarding the treatment 
of the time/energy dependence of the key quantities, is 
in progress to reduce the prefactor of the scaling of our 
method. 
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VII. SUMMARY 

We have presented a detailed account of a method for 
calculating the electron self-energy and related quanti- 
ties from ab-initio many-body perturbation theory within 
the GW approximation. The method is based on repre- 
senting the basic quantities (Green's function, dielectric 
response function, dynamically screened Coulomb inter- 
action and self-energy) on a real-space grid and on the 
imaginary time axis. In those intermediate steps of the 
calculation where it is computationally more efficient to 
work in reciprocal space and imaginary energy we change 
to the latter representation by means of fast Fourier 
transforms. Working on the imaginary time/energy axis 
considerably facilitates the numerical treatment. The 
matrix elements of the self-energy on the real-energy axis 
are obtained by an analytic continuation procedure which 
was shown to be accurate and stable. We have demon- 
strated the accuracy of the method by calculating quasi- 
particle excitation energies at high-symmetry points for 
the prototype semiconductor Si. The computational ef- 
fort of the method scales quadratically with unit cell size. 
This was shown explicitly by calculations for Si super- 
cells. The method allows the extension of ab-initio work 
beyond the calculation of quasiparticle energies and its 
application to materials with larger basis sets or larger 
unit cells than were previously feasible. Calculating the 
full self-energy gives access to quantities like the charge 
density p^ , spectral function p^ , momentum distribu- 
tion and total energy at the GW level. It also opens the 
possibility of calculating the self-energy self-consistently 
and to provide the GW Green's function (after solving 
the Dyson equation) - essential prerequisites for going be- 
yond the GW approximation, e. g., by iterating Hedin's 
equations. This remains to be investigated in future. 
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FIG. 1. The grids and cells used in the calculation in (a) 
real and (b) reciprocal space. These are shown here schemat- 
ically in two dimensions for the case of a 2 x 2 grid of k-points 
(corresponding to an interaction cell of 2 x 2 unit cells), and a 
3x3 real-space grid in the unit cell (corresponding to a 3 x 3 
grid of reciprocal lattice vectors in reciprocal space) . The bare 
Coulomb interaction 1/ |r — r'| is non-periodic on the interac- 
tion cell, and its amplitude at a real-space-grid point is taken 
to be that at the corresponding point in the Wigner-Seitz cell 
around r of the IC lattice. All other quantities, in both real 
and reciprocal space, are periodic. In these cases the choice 
of the primitive (shown) or Wigner-Seitz cell is a matter of 
computational convenience. 



FIG. 2. Fitting of the matrix elements of the correlation 
self-energy (Ec(iw)) of silicon as a function of imaginary en- 
ergy. The real (top panel) and imaginary (lower panel) parts 
are both reproduced extremely accurately by the fits (which 
on this scale are indistinguishable from the calculated values). 
The bands shown are band 1 at the F point (solid line for fit- 
ted and calculated function), bands 2-4 at F (dotted line), 
bands 3-4 at X (dashed) and bands 5-6 at X (dot-dashed). 



FIG. 3. Convergence of a matrix element (degenerate 
bands 2-4 at F) of the correlation self-energy (Ec(ici^)) of sili- 
con with respect to Umax with a fixed energy grid spacing of 
Auj= 0.16 Har. The top panel shows the calculated self-energy 
on the imaginary axis with the analytically continued depen- 
dence on real energy shown in the lower panel. The lines 
correspond to LUmax= 5 Har (dotted), 10 Har (dashed) and 
20 Har (solid). Changes in {Ec(a;)) are directly related to 
changes in the calculated (Ec(iaj)), rather than any instabili- 
ties in the analytic continuation technique. 



FIG. 4. Convergence of (Ec(jci^)) with respect to energy 
grid spacing Alj with fixed Umax= 10 Har. The lines corre- 
spond to /S.u>^ 0.24 Har (dotted), 0.16 Har (dashed) and 0.08 
Har (solid). As in the previous Figure the analytically contin- 
ued matrix element converges well, indicating the stability of 
the fitting procedure. The matrix elements on the real-energy 
axis (which are obtained in form of fits to a model function, 
see text) are plotted with a fixed grid spacing of 0.04 Har in 
order to facilitate comparison between the three curves. 



FIG. 5. Real part of a matrix element (degenerate bands 
2-4 at F) of the correlation self-energy (Ec(ti;)) on the real en- 
ergy axis, computed with Umax = 5 Har (dot-dashed line), 10 
Har (dashed line), and 20 Har (long-dashed line) with fitting 
{Ec(iaj)) on the imaginary energy axis using the energy range 
{~Umax,+Umax) and With Umax = 10 Har (dotted fine) and 
20 Har (solid line) using a fixed energy range of (—5 Har, -1-5 
Har) for the fitting. 



FIG. 6. Scaling of the CPU time on a Digital Alpha 
500/500 workstation with respect to unit cell size for Sin 
(n = 2,4,8). Besides the total CPU time (o) the times for 
three major parts of the computation are given: calculation 
of (i) the Green's function (□), (ii) the dynamically screened 
Coulomb interaction (O), and (iii) the matrix elements of the 
self-energy including a recomputation of the Green's function 
(A), see text. Note the quadratic spacing of the abscissa. 



TABLE I. Convergence of quasiparticle energies for bulk 
Si with respect to cutoff and grid parameters in the present 
method. Shown are the parameters needed for a convergence 
of 50 meV and 20 meV, respectively. 



parameter 



50 meV 



20 meV 



plane-wave cutoff (in Ry) 13.5 16 

band cutoff (in Ry) 8 10 

size of k grid 4x4x4 5x5x5 
spacing /S.t of time grid (in a.u.) 0.3 0.15 
range tmax of time grid (in a.u.) 13. 20. 
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TABLE II. Calculated quasiparticle energies at points of 
high symmetry for Si (in eV). The valence band maximum has 
been set to zero. The calculation was done with a plane- wave 
cutoff of 16 Ry, a 5 X 5 X 5 k-point grid, a band cutoff of 
10 Ry and a time grid with At = 0.3 a.u. and Tmax = 20 
a.u. For comparison the eigenvalues obtained in the LDA 
calculation providing the input for the GW calculation have 
been included. GW 1 refers to a GW calculation using the 
LDA Green function whereas in GW 2 the Green function 
has been updated by replacing the LDA eigenvalues by the 
quasiparticle energies obtained in GW 1. In the last column 
some experimental data are given (see text). 



band 



Ti,, 

r' 

r' 

J- 15c 

X\c 
X4c 

T' 

^2v 

Llv 

Lie 

L2c 



LDA 



-11.89 
0.00 
2.58 
3.28 

-7.78 
-2.82 
0.61 
10.11 

-9.57 
-6.96 
-1.17 
1.46 
3.33 
7.71 

0.49 



GW 1 



GW 2 



-11.57 
0.00 
3.24 
3.94 

-7.67 
-2.80 
1.34 
10.54 

-9.39 
-6.86 
-1.17 
2.14 
4.05 
8.29 

1.20 



-11.58 
0.00 
3.32 
4.02 

-7.68 
-2.81 
1.42 
10.63 

-9.40 
-6.88 
-1.17 
2.22 
4.14 
8.39 

1.28 



Experiment 



-12.5 ±0.6'' 
0.00 



3.40=^, 
4.23=" 



3.05'= 
, 4.1'' 



-2.9^ -3.3 ±0.2=* 
1.25*^ 

-9.3±0.4'- 
-6.7 ±0.2^ 
-1.2 ±0.2=^, -1.5° 
2.1* , 2.4 ±0.15*^ 
4.15 ±0.1*5 

1.17=^ 



TABLE IV. Comparison of quasiparticle energies (in eV) 
at high symmetry points for bulk Si with the corresponding 
values obtained from calculations for Si4 and Sig supercells 
(see text). All energies refer to the respective valence band 
top. Shown are the band energies for all k-points mapping 
onto r, X, and L, respectively, when the unit cell size is 
increased, where those k-points appear in the calculation. 







bulk Si 


Si4 




Sig 




r 


-11.55 


-7.63 -10.53 


-11.53 -7.61 


-11.55 


-7.63 


-10.53 




0.00 


-2.78 -3.41 


0.00 -2.78 


0.00 


-2.78 


-3.41 




3.23 


1.35 -1.81 


3.22 1.37 


3.22 


1.37 


-1.81 




3.97 


10.52 1.77 


3.96 10.45 


3.96 


10.43 


1.78 






3.82 








3.82 






6.35 








6.30 


L 


-9.35 


-9.62 


-9.34 


-9.36 


-9.63 






-6.83 


-5.35 


-6.81 


-6.83 


-5.35 






-1.16 


-3.62 


-1.16 


-1.16 


-3.62 






2.17 


-1.27 


2.18 


2.18 


-1.27 






4.04 


3.13 


4.03 


4.02 


3.12 






8.28 


5.46 


8.26 


8.25 


5.42 








5.66 






5.63 








6.75 






6.70 




X 


-7.63 


-7.46 


-7.61 


-7.64 


-7.47 






-2.78 


-3.73 


-2.78 


-2.79 


3.73 






1.35 


4.84 


1.37 


1.36 


4.82 






10.52 


5.67 


10.44 


10.42 


5.63 





"Ref. 
''Ref. 
'^Ref. 
''Ref. 
"Ref. 



32|. 

M. 



'Ref. |36 
=Ref. 



TABLE III. Grid parameters for the Si„ (n = 2, 4, 8) cal- 
culations (see text). 



parameter 


Si2 

(bulk Si) 


Si4 


Sig 


Niuc 


55 


230 


455 


Nuc 


9x9x9 


9 X 9 X 18 


9 X 9 X 36 


Nk 


4x4x4 


4x4x2 


4x4x1 


Nic = NucNk 


36 X 36 X 36 


36 X 36 X 36 


36 X 36 X 36 




63 


63 


63 
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(a) Real Space 




(b) Reciprocal Space 
Origin 
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